Trapping oscillations, discrete particle effects and kinetic theory of collisionless plasma 



O 

o 

(N 



in ■ 
O ■ 

Oh- 



in 

(N 
O 
cn 
o 



CZ3 

o 

^ : 

Oh. 



X 



F. Dovcil'^i, M-C. Firpo'^lil, Y. Elskens^lll, D. Guyomarc'h'^, M. Polcni'' and P. Bertrand'S 

Eqmpe turbulence plasma, Physique des interactions ioniques et moleculaires, 
Unite 6633 CNRS-Universite de Provence, 
case 321, Centre de Saint- Jerome, F-13397 Marseille cedex 20 
'' Laboratoire de physique des milieux ionises et applications, 
Unite 7040 CNRS-Universite H. Poincare, Nancy I, 
BP 239, F-545O6 Vandceuvre cedex, France 
(preprint TP99.il) 



Effects induced by the finite number A'^ of particles on tlie 
evolution of a monochromatic electrostatic perturbation in a 
collisionless plasma are investigated. For growth as well as 
damping of a single wave, discrete particle numerical simula- 
tions show a A'^-dependent long time behavior which differs 
from the numerical errors incurred by vlasovian approaches 
and follows from the pulsating separatrix crossing dynamics 
of individual particles. 
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I. INTRODUCTION 

It is tempting to expect that kinetic equations and 
their numerical simulation provide a fair description of 
the time evolution of systems with long-range or 'global' 
interactions. A typical, fundamental example is offered 
by wave-particle interactions, which play a central role 
in plasmas. In this Letter we test this opinion explicitly. 

Collisionless plasma dynamics is dominated by collec- 
tive processes. Langmuir waves and their familiar Lan- 
dau damping and growth Q are a good example of these 
processes, with many applications, e.g. plasma heating 
in fusion devices and laser-plasma interactions. For sim- 
plicity we focus on the one-dimensional electrostatic case, 
traditionally described by the (kinetic) coupled set of 
Vlasov-Poisson equations The current debate on 

the long-time evolution of this system shows that further 
insight in this fundamental process is still needed 

The driving process (induced by the binary Coulomb 
interaction between particles) is the interaction of the 
electrostatic waves in the plasma with the particles at 
nearly resonant velocities, which one analyses canonically 
by partitioning the plasma in bulk and tail particles. 



Langmuir modes are the collective oscillations of bulk 
particles, with slowly varying complex amplitudes in an 
envelope representation ; their interaction with tail par- 
ticles is described by a self-consistent set of hamiltonian 
equations pi. These equations already provided an effi- 
cient basis |6|] for investigating the cold beam plasma in- 
stability and exploring the nonlinear regime of the bump- 
on-tail instability Q . Analytically, they were used to give 
an intuitive and rigorous derivation of spontaneous emis- 
sion and Landau damping of Langmuir waves Be- 
sides, as it eliminates the rapid plasma oscillation scale 
ojp^, this self-consistent model offers a genuine tool to 
investigate long-time dynamics. 

As we follow the motion of each particle, we can also 
address the influence of the finite number of particles in 
the long run. This question is discarded in the kinetic 
Vlasov-Poisson description, for which the finite-iV cor- 
rection is the Balescu-Lenard equation ||^ formally de- 
rived from the accumulation of weak binary collisions, 
with small change of particle momenta. It implies a dif- 
fusion of momenta, driving the plasma towards equilib- 
rium. However, when wave-particle coupling is domi- 
nant, the Balescu-Lenard equation is not a straightfor- 
ward approach to finite- effects on the wave evolution. 

Here we investigate direct finite- A^ effects on the self- 
consistent wave-particle dynamics. It is proved | pO| that 
the kinetic limit N 00 commutes with evolution over 
arbitrary times. As one might argue that finite A^ be 
analogous to numerical discretisation in solving kinetic 
equations, we also integrate the kinetic system with a 
'noise-free' semi-lagrangian solver In this Letter we 
compare finite grid effects of the kinetic solver and gran- 
ular aspects of the A^-particles system, whose evolution 
is computed with a symplectic scheme . 

We discuss the case of one wave interacting with the 
particles. Though a broad spectrum of unstable waves is 
generally excited when tail particles form a warm beam, 
the single-wave situation can be realized experimentally 
]l2t and allows to leave aside the difficult problem of 
mode coupling mediated by resonant particles |13[. 
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II. SELF-CONSISTENT WAVE-PARTICLE 
MODEL AND KINETIC MODEL 

Consider a one-dimensional electrostatic potential per- 
turbation ^{z, t) = [<j)k{T) exp i{kz — uj^t) + c.c] (where 
c.c. denotes complex conjugate), with complex enve- 
lope (pki in a plasma of length L with periodic boundary 
conditions (and neutralizing background). Wavenum- 
ber k and frequency LOk satisfy a dispersion relation 
e{k,ujk) = 0. The density of N (quasi-)resonant elec- 
trons is (t(z,t) — inL/N)J2i=iS{z ~ ^^/(t)), where n is 
the electron number density and zi is the position at time 
T of electron labeled I (with charge e and mass m) . Non- 
resonant electrons contribute only through the dielectric 
function e, so that ipk and the zj's obey coupled equations 



d(t>k/dT 



N 



eok^N{de/dujk) j-[ 



exp[-i/cz; + iuJkT] (1) 



d^zi/dr^ = {iek/m)(f>kex.]i[ikzi — iw^r] -fee. (2) 

where eg is the vacuum dielectric constant. With = 
ne^/[meo(9e/9cjfc)] [p|, i = ar, ' = d/dt, xi = kzi - LUkT 
and V = (efc^0fe)/(a^m), this system defines the self- 
consistent dynamics (with iV -f 1 degrees of freedom) 



N 



-IXl) 



xi = zT^exp(zx/) — iV* exp(— ix;) 



(3) 
(4) 



for the coupled evolution of electrons and wave in di- 
mensionless form. This system derives from hamil- 
tonian i7^x,p,C,C) = E/li(P?/2 - N-^/\e'-^ - 
N ^/^C*e "'), where a star means a complex conju- 
gate and C = N'^^'^V . An efficient symplectic integration 
scheme is used to study this hamiltonian numerically . 

The system (|)-(|) is invariant under two continuous 
groups of symmetries. Invariance under time translations 
implies the conservation of the energy V\ = H . The phase 
9 of Q = |C|e"*^ plays the role of a position for the wave, 
and system (ph-(y) is also invariant under translations 



^ 6 -\- a, Xi = xi 



This translation invariance leads 



to the conservation of momentum P = pi + where 
the contribution from the wave is analogous to the Poynt- 
ing vector of electromagnetic waves (which is quadratic 
in the electromagnetic fields) [|l6| . Conservation of these 
invariants constrains the evolution of our system, and we 
checked that the numerical integration preserves them. 

In the kinetic limit iV — > oo, electrons are distributed 
with a density f{x,p,t), and system (||)-(||) yields the 
Vlasov-wave system 

V = i [ e-"f{x,p,t)dxdp (5) 



For initial data approaching a smooth function / as 
N CO, the solutions of (^-(^ converge to those of 
the Vlasov-wave system over any finite time interval [ p^ . 
This kinetic model is integrated numerically by a semi- 
lagrangian solver, covering {x,p) space with a rectangular 
mesh : the function / (interpolated by cubic splines) is 
transported along the characteristic lines of the kinetic 
equation, i.e. along trajectories of the original particles 

Let us first study linear instabilities. One solution of 
(^)-(^) corresponds to vanishing field Vq = 0, with par- 
ticles evenly distributed on a finite set of beams with 
given velocities. Small perturbations of this solution have 
SV = SVoe''^, with rate 7 solving [|| 



7 = 7r 4- 271 = iN 



N 

-1 \ ^ 

1=1 



{l + ipiY 



(J) 



For a monokinetic beam with velocity U , (0) reads 7(7 -f 
iUY = i ; flic most unstable solution occurs for [7 = 
(with 7 = (V3-fi)/2). For a warm beam with smooth 
initial distribution f{p) (normalized to J fdp = 1), the 
continuous limit of (^ yields 7 = i Jij + ip)^^ f{p)dp. 
For a sufficiently broad distribution (|/'(0)| ^ 1), we 
obtain |7i|7r — 7i'7r/'(— 7i), where /' = df/dp, and 
7i ~ 7r7r/"(0) for |/"(0)| < tt-^. Except for the triv- 
ial solution 7i. = 0, other solutions can only exist for 
a positive slope /'(O). Then the perturbation is unsta- 
ble as the evolution of SV is controlled by the eigen- 
value 7 with positive real part, i.e. with growth rate 
7r « 7l = 7i'/'(0) > 0. Negative slope leads to the lin- 
ear Landau damping paradox : the observed decay rate 
7l — 7r/'(0) < is not associated to genuine eigenvalues, 
but to phase mixing of eigenmodes [p| jr^ , p^ , as a direct 
consequence of the hamiltonian nature of the dynamics. 

Now, this linear analysis generally fails to give the large 
time behavior. This is obvious for the unstable case as 
non-linear effects are no longer negligible when the wave 
intensity grows so that the trapping frequency tOhit) = 
y^2\V{t)\ becomes of the order of the linear rate 7^. 

We used the monokinetic case as a testbed |18| ; |l9| . 
Finite-A^ simulations show that the unstable solution 
grows as predicted and saturates to a limit-cycle-like be- 
havior where the trapping frequency ujh{t) oscillates be- 
tween 1.271- and 271-. In this regime, some of the ini- 
tially monokinetic particles have been scattered rather 
uniformly over the chaotic domain, in and around the 
pulsating resonance, while others form a trapped bunch 
inside this resonance (away from the separatrix) p9| . 
This dynamics is quite well described by effective hamil- 
tonians with few degrees of freedom 
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dtf + pdj + {iVe'^ - iV*e-ndpf = (6) 



In this Letter, we discuss the large time behavior of the 
warm beam case, with f'{po) ^ at the wave nominal 
velocity po = 0. Fig. ^ displays three distribution func- 
tions (in dimensionless form) with similar velocity width : 
(i) & function (CD) giving the same decay rate for all 
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phase velocities, (ii) a function (CG) giving a constant 
growth rate for all phase velocities 0, (Hi) a truncated 
Lorentzian (TL) with positive slope /'(O) > 0. 

III. DAMPING CASE 

For the damping case, the linear description introduces 
time secularities which ultimately may break linear the- 
ory down : the ultimate evolution is intrinsically nonlin- 
ear, not only if the initial field amplitude is large, as in 
O 'Neil's seminal picture ||^, but also if one considers the 
evolution over time scales of the order of the trapping 
time (which is large for small initial wave amplitude). 
The question of the plasma wave long-time fate is thus 
far from trivial Though some simulations infer 
that nonlinear waves eventually approach a Bernstein- 
Greene-Kruskal steady state instead of Landau van- 
ishing field, the answer should rather strongly depend on 
initial conditions. Our TV-particle, 1-wave system is the 
simplest model to test these ideas. 

A thermodynamical analysis predicts that, for a 
warm beam and small enough initial wave amplitude, 
ujh ^ N~^/'^ at equilibrium in the limit ^ oo. Fig. ^ 
shows the evolution of a small amplitude wave launched 
in the beam. The A-particle system (line N) and the 
kinetic system (line V) initially damp the wave exponen- 
tially as predicted by perturbation theory |^ , for a time 
of the order of |7l|~^- 

After that phase-mixing time, trapping induces non- 
linear evolution and both systems evolve differently. For 
the A-particle system, the wave grows to a thermal level 
that scales as A~^/^, corresponding to a balance be- 
tween damping and spontaneous emission p|,p^. For 
the kinetic system, initial Landau damping is followed 
by slowly damped trapping oscillations around a mean 
value which also decays to zero, at a rate decreasing for 
refined mesh size. Fig. ^ reveals that finite- A^ and kinetic 
behaviors can considerably diverge as spontaneous emis- 
sion is taken into account. The time t^t after which the 
finite-A effects force this divergence is found to diverge 
as A CXI. 

IV. UNSTABLE CASE 

Now consider an unstable warm beam (/'(O) > 0). 
Line Nl (resp. N2) of Fig. ^ displays ln(a;b(t)/7r) versus 
time for (|)-(|) with a CG distribution with N = 128000 
(resp. 512000) and 7^ = 0.08. Line VI (resp. V2) shows 
ln(a;b(t)/7r) versus 7ri for the kinetic system and the 
same initial distribution with a 32 x 128 (resp. 256 x 1024) 
grid in {x,p) space. All four lines exhibit the same initial 
exponential growth of linear theory with less than 1% er- 
ror on the growth rate. Saturation occurs for Wb/7r ~ 3.1 
1^. Lines Nl and VI do not superpose beyond the first 



trapping oscillation after saturation. Note that, in our 
system, oscillating saturation does not excite sideband 
Langmuir waves as our hamiltonian incorporates only a 
single wave, not a spectrum. 

After the first trapping oscillation, kinetic simulations 
exhibit a second growth at a rate controlled by mesh size. 
Line V2 suggests that a kinetic approach would predict 
a level close to the trapping saturation level on a time 
scale awarded by reasonable integration time. This level 
is fairly below the equilibrium Vth predicted by a gibbsian 
approach ; such pathological relaxation properties 
in the A — > 00 limit seem common to mean-field long- 
range models p3|] . Both kinetic simulations also exhibit 
a strong damping of trapping oscillations, which disap- 
pear after a few oscillations, whereas finite-A simulations 
show persistent trapping oscillations. 

One could expect that finite-A effects would mainly 
damp these oscillations, so that the wave amplitude 
reaches a plateau. Actually, we observe persistent os- 
cillations for all A, and the wave amplitude slowly grows 
further, whereas the velocity distribution function fiat- 
tens over wider intervals of velocity. 

This spreading of particles is due to separatrix cross- 
ings, i.e. successive trapping and detrapping by the wave 
jrof . Indeed, when the wave amplitude grows (during 
its pulsation), it captures particles with nearby veloc- 
ity, i.e. with a relative velocity Awjn « i^/SlV^ ; the 
trapped particles start bouncing in the wave potential 
well. When the wave amplitude decreases, particles are 
released, but if they experienced only half a bouncing 
period, they are released with a relative velocity (with 
respect to the wave) opposite to their initial one, i.e. 
Auout ~ — Atijn. Now notice that a particle which has 
just been trapped would oscillate at a longer period than 
the nominal bouncing period (namely the one deep in 
the potential). Moreover, if the recently trapped particle 
had just adiabatic motion in the well, it would have to 
recross the separatrix when the resonance would enclose 
the same area as at its trapping |Q. Thus one expects 
the particle to be unable to complete a full bounce, and 
the fraction of particles for which Awout ~ —^V'm is sig- 
nificant. 

During this particle spreading process in (x,p) space, 
the wave pulsation is maintained by the bunch of parti- 
cles which were initially trapped, and are deep enough in 
the potential well to remain trapped over a whole bounc- 
ing period. These particles form a macroparticle, as is 
best seen in the case of a cold beam Note that, over 
long times, the macroparticle must slowly spread in the 
wave resonance, following two processes. One acts if the 
trapped particle motion is regular : the trapped motions 
are anisochronous, i.e. have different periods (only the 
harmonic oscillator has isochronous oscillations). The 
other one works if the motion is chaotic : nearby trajec- 
tories diverge due to chaos. Both processes contribute to 
the smoothing of the particle distribution for long times. 
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but over much longer times than those over which we 
follow the system evolution and observe the wave modu- 
lation. 

This second growth after the first trapping saturation 
depends on the shape of the initial distribution function. 
In Fig. ||(b), line N2 is the same as in Fig. ^(a), com- 
puted over a longer duration, and line N3 corresponds to 
N = 64000 with the TL distribution of Fig. Although 
N3 corresponds to 8 times fewer particles than N2, the 
final level reached at the end of the simulation is lower. 
In the second growth, particles are transported further 
in velocity, so that the plateau in f{p) broadens with 
time. As the wave grows, it can trap particles with ini- 
tial velocity further away from its phase velocity. Since 
the TL distribution reaches its maximum at w « 1.06 and 
decreases significantly beyond this velocity (while CG is 
still growing for larger v), fewer particles (with TL than 
with CG) can give momentum to the wave when being 
trapped (P is conserved) ; hence the second growth is 
slower for the TL distribution. 

We followed the evolution of the wave amplitude for N3 
up to 7i.t = 1750 : starting from the first trapping sat- 
uration level (0.4Vth), fluctuations persist with a growth 
rate that slowly decreases as we reach 0.78I4h at the end 
of the computation. Line N4 of Fig. || corresponds to the 
TL distribution with 2048000 particles and shows persis- 
tent oscillations with approximately the same amplitude 
as for TV ^ 64000. 

V. CONCLUSION 

These observations clearly indicate that the kinetic 
models are an idealization and do not contain all the 
intricate behavior of a discrete particles system. Now, 
we must also admit that the kinetic simulation schemes 
do not exactly reproduce the analytic implications of the 
kinetic equation. It is then legitimate to ask whether the 
numerical implementation of the kinetic equations repro- 
duce the difference between the finite- dynamics and 
the kinetic theory. 

A basic property of the coUisionless kinetic equation is 
that it transports the distribution function f{x,p) along 
the particle trajectories (or characteristic lines in (x,p) 
space). As long as the kinetic calculation of / is accu- 
rate, one expects the kinetic scheme to follow closely the 
A^-particle dynamics too. However, the kinetic scheme 
is bound to depart from the analytic predictions of the 
kinetic equation, because the (chaotic or anisochronous) 
separation of particle trajectories implies that constant- 
/ contours eventually evolve into complex, interleaved 
shapes. This filamentation is smoothed by numerical par- 
tial differential equation integrators, while A^-body dy- 
namics follows the particles more realistically, sustaining 
the trapping oscillations. Hence both types of dynamics 



will depart from each other when filamentation reaches 
scales below the semi-lagrangian kinetic code grid mesh. 

The onset of filamentation is easily evidenced in kinetic 
simulations. Indeed, whereas the kinetic equation analyt- 
ically preserves the 2-entropy /(I — f)fdxdp, numerical 
schemes increase entropy significantly when constant-/ 
contours form filaments in (a;,p)-space As this is 

also the time at which trapping oscillations are found to 
damp in our simulations, it appears that vlasovian sim- 
ulations must be considered with caution from that time 
on - and it turns out that it is also the time from which 
the second growth starts. 

In summary, discussing the basic propagation of a sin- 
gle electrostatic wave in a warm plasma, we presented 
finite-A^ effects which do not merely result from nu- 
merical errors and elude a kinetic simulation approach. 
Their understanding depends crucially on the dynamics 
in phase space. The sensitive dependence of microscopic 
evolution to the fine structure of the initial particle dis- 
tribution in phase space [l^ implies that the interplay 
between limits t —^ oo and N ^ oo requires some cau- 
tion. Somewhat paradoxically, refining the grid for the 
Vlasov simulations does not solve this problem. 

The driving process in the system evolution is sepa- 
ratrix crossing, which requires a geometric approach to 
the system dynamics. Further work in this direction p^ ] 
will also shed new light on the foundations of common 
approximations, such as replacing original dynamics (|l|)- 
(^) by coupled stochastic equations, in which particles 
undergo noisy transport. 
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FIG. 2. Time evolution of ln(a'b(t)/|7L|) for a CD velocity 
distribution and initial wave amplitude below thermal level : 

(N) A'^-particles system with A'^ — 32000, (V) kinetic scheme 
with 32 X 512 {x,p) grid. Inset : short-time evolution. 
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FIG. 3. Time evolution of In (t<;b(t)/7r). (a) CG initial dis- 
tribution : kinetic scheme with (VI) 32 x 128, (V2) 256 x 1024 
(.T,p) grid ; A^-particlcs system with (Nl) N = 128000, (N2) 
N = 512000 ; (b) Comparison of CG (N2) with TL initial 
distribution for (N3) N = 64000, (N4) N = 2048000. 
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